library(sf)
library(readxl)
library(dplyr)
library(plyr)
library(ggplot2)
library(afrihealthsites)
library(ggpubr)
library(afriadmin)
library(tmap)
library(cowplot)
library(forcats)
library(colorspace)
library(readr)
library(sunburstR)
library(htmltools)
library(htmlwidgets)

# Install Malawi MFL
malawi_MFL = read_excel("~/malawi-health-facilities-1/MHFR_Facilities 1.xlsx")

# Convert to sf

## omit NA's
new_malawi_MFL = na.omit(malawi_MFL)

## check for NA 
any(is.na(new_malawi_MFL))
## [1] FALSE
## transform geometry columns into numeric 
new_malawi_MFL = transform(new_malawi_MFL, LATITUDE = as.numeric(LATITUDE), 
                                           LONGITUDE = as.numeric(LONGITUDE))

## convert to sf object
malawi_facilities_MFL = st_as_sf(new_malawi_MFL, coords = c("LONGITUDE", "LATITUDE"), dim = "XY")

malawi_facilities_MFL = st_set_crs(malawi_facilities_MFL, 4326) ## set CRS

head(malawi_facilities_MFL)
## Simple feature collection with 6 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 33.74129 ymin: -15.84 xmax: 35.09 ymax: -13.79742
## geographic CRS: WGS 84
##       CODE                                               NAME  COMMON.NAME
## 1 MC010002                               A + A private clinic          A+A
## 2 BT240003                                       A-C Opticals A.C Opticals
## 3 BT240005                                Akwezeke PVT Clinic Akwezeke Pvt
## 4 BT240006                                  AB Medical Clinic        Abowa
## 5 LL040007                                 ABC Comm. Hospital   ABC Clinic
## 6 LL040010 Achikondi Women Community Friendly Services Clinic    Achikondi
##                                       OWNERSHIP       TYPE     STATUS
## 1                                       Private     Clinic Functional
## 2                                       Private     Clinic Functional
## 3                                       Private     Clinic Functional
## 4                                       Private     Clinic Functional
## 5 Christian Health Association of Malawi (CHAM)   Hospital Functional
## 6                                       Private Dispensary Functional
##                 ZONE DISTRICT DATE.OPENED                   geometry
## 1 Centrals West Zone  Mchinji  Jan 1st 75 POINT (33.88563 -13.79742)
## 2    South East Zone Blantyre  Jan 1st 75        POINT (35.03 -15.8)
## 3    South East Zone Blantyre  Jan 1st 75       POINT (35.09 -15.84)
## 4    South East Zone Blantyre  Jan 1st 75       POINT (35.09 -15.84)
## 5 Centrals West Zone Lilongwe  Jan 1st 75 POINT (33.74129 -13.96816)
## 6 Centrals West Zone Lilongwe  Jan 1st 75  POINT (33.7793 -13.95473)

Overview:

  1. Malawi MFL
# Re-order the facility types
malawi_facilities_MFL$TYPE = as.factor(malawi_facilities_MFL$TYPE)

malawi_facilities_MFL$TYPE = factor(malawi_facilities_MFL$TYPE, levels = c("Central Hospital", "District Hospital", "Hospital", "Health Centre", "Clinic", "Health Post", "Dispensary", "Private", "Unclassified"))

# Number of each type + ownership
facility_types_MFL = as.data.frame(table(malawi_facilities_MFL$TYPE, malawi_facilities_MFL$OWNERSHIP))
head(facility_types_MFL)
##                Var1            Var2 Freq
## 1  Central Hospital Aquaid Lifeline    0
## 2 District Hospital Aquaid Lifeline    0
## 3          Hospital Aquaid Lifeline    0
## 4     Health Centre Aquaid Lifeline    0
## 5            Clinic Aquaid Lifeline    0
## 6       Health Post Aquaid Lifeline    0
levels(facility_types_MFL$Var2)
## [1] "Aquaid Lifeline"                              
## [2] "Christian Health Association of Malawi (CHAM)"
## [3] "Government"                                   
## [4] "Mission/Faith-based (other than CHAM)"        
## [5] "Non-Government"                               
## [6] "Other"                                        
## [7] "Parastatal"                                   
## [8] "Private"
# Only showing government, private and CHAM ownership
facility_types_MFL$new_Var2 = revalue(facility_types_MFL$Var2, c("Aquaid Lifeline"="Other", "Non-Government"="Other", "Parastatal"="Other", "Mission/Faith-based (other than CHAM)"="Other"))
facility_types_MFL$new_Var2 = factor(facility_types_MFL$new_Var2, levels = c("Government", "Private", "Christian Health Association of Malawi (CHAM)", "Other"))


## bar plot of no. of facility types 
plot_facility_types_MFL = ggplot(facility_types_MFL, aes(x=Var1, y=Freq, fill=new_Var2)) + geom_bar(position = "stack", stat = "identity")

plot_facility_types_MFL = plot_facility_types_MFL + labs(x = "Facility types", y = "Frequency", fill="Ownership") + scale_fill_brewer(palette = "Set2") + coord_flip() + theme_minimal() + ggtitle("MFL") + theme(axis.title.y = element_blank(), legend.title = element_text(size = 8, face = "bold"), legend.text = element_text(size = 8), plot.title = element_text(face = "bold", hjust = 0), legend.position = "bottom")
plot_facility_types_MFL

par(mar=c(11,4,4,4))
# Re-order ownership
malawi_facilities_MFL$OWNERSHIP = as.factor(malawi_facilities_MFL$OWNERSHIP)

malawi_facilities_MFL$OWNERSHIP = factor(malawi_facilities_MFL$OWNERSHIP, levels = c("Government", "Private", "Christian Health Association of Malawi (CHAM)", "Non-Government", "Mission/Faith-based (other than CHAM)", "Other", "Parastatal", "Aquaid Lifeline"))

# Number of each type of ownership
ownership_MFL = as.data.frame(table(malawi_facilities_MFL$OWNERSHIP))

## bar plot of ownership
plot_ownership_MFL = ggplot(ownership_MFL, aes(x=forcats::fct_relabel(Var1,stringr::str_wrap,width = 16), y=Freq)) + geom_bar(stat = "identity", fill="slategray")

plot_ownership_MFL = plot_ownership_MFL + labs(x="Ownership", y = "Frequency") + coord_flip() + theme_minimal() + ggtitle("MFL") + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold"))
plot_ownership_MFL

  1. WHO-KWTRP
# Malawi WHO data.frame
malawi_WHO <- afrihealthsites("malawi", datasource='who', plot=FALSE, returnclass='dataframe')

head(malawi_WHO)
## # A tibble: 6 x 10
##   Country Admin1 `Facility name` `Facility type` Ownership   Lat  Long
##   <chr>   <chr>  <chr>           <chr>           <chr>     <dbl> <dbl>
## 1 Malawi  Centr… 80 Block Clinic Clinic          MoH       -12.9  33.4
## 2 Malawi  Centr… ABC Community … Clinic          FBO       -14.0  33.7
## 3 Malawi  Centr… Adventist Heal… Health Centre   FBO       -14.0  33.8
## 4 Malawi  Centr… Alinafe Commun… Community Hosp… FBO       -13.4  34.2
## 5 Malawi  Centr… Area 18 Health… Health Centre   MoH       -13.9  33.8
## 6 Malawi  Centr… Area 25 Health… Health Centre   MoH       -13.9  33.8
## # … with 3 more variables: `LL source` <chr>, iso3c <chr>,
## #   facility_type_9 <chr>
## convert malawi_WHO to sf object

class(malawi_WHO)
## [1] "tbl_df"     "tbl"        "data.frame"
any(is.na(malawi_WHO))
## [1] TRUE
new_malawi_WHO = na.omit(malawi_WHO) ## omit NA

sf_malawi_WHO = st_as_sf(new_malawi_WHO, coords = c("Long", "Lat"), dim = "XY")
sf_malawi_WHO = st_set_crs(sf_malawi_WHO, 4326)
# Re order facility types and ownership 
sf_malawi_WHO$`Facility type` = as.factor(sf_malawi_WHO$`Facility type`)
sf_malawi_WHO$`Facility type` = factor(sf_malawi_WHO$`Facility type`, levels = c("Central Hospital", "District Hospital", "Mission Hospital", "Rural Hospital", "Community Hospital", "Health Centre", "Clinic", "Health Post/Dispensary"))

sf_malawi_WHO$Ownership = as.factor(sf_malawi_WHO$Ownership)
sf_malawi_WHO$Ownership = factor(sf_malawi_WHO$Ownership, levels = c("MoH", "FBO", "Local authority", "NGO"))

# No. of original facility types + ownership
facility_types_WHO = as.data.frame(table(sf_malawi_WHO$`Facility type`, sf_malawi_WHO$Ownership))
head(facility_types_WHO)
##                 Var1 Var2 Freq
## 1   Central Hospital  MoH    4
## 2  District Hospital  MoH   24
## 3   Mission Hospital  MoH    0
## 4     Rural Hospital  MoH   17
## 5 Community Hospital  MoH    0
## 6      Health Centre  MoH  330
## bar plot of original facility types 
plot_facility_types_WHO = ggplot(facility_types_WHO, aes(x=Var1, y=Freq, fill=Var2)) + geom_bar(position = "stack", stat = "identity")
plot_facility_types_WHO = plot_facility_types_WHO + labs(x = "Facility types", y = "Frequency", fill="Ownership") + scale_fill_brewer(palette = "Set2") + coord_flip() + theme_minimal() + ggtitle("WHO") + theme(axis.title.y = element_blank(), legend.title = element_text(size = 8, face = "bold"), legend.text = element_text(size = 8), plot.title = element_text(face = "bold", hjust = 0), legend.position = "bottom") + expand_limits(y=c(0,600))
plot_facility_types_WHO

# Re order
sf_malawi_WHO$facility_type_9 = as.factor(sf_malawi_WHO$facility_type_9)
sf_malawi_WHO$facility_type_9 = factor(sf_malawi_WHO$facility_type_9, levels = c("Hospital", "Health Centre", "Health Clinic", "Health Post", "Community Health Unit"))

# No. of reclassified facility types
RC_facility_types_WHO = as.data.frame(table(sf_malawi_WHO$facility_type_9))
RC_facility_types_WHO
##                    Var1 Freq
## 1              Hospital   80
## 2         Health Centre  451
## 3         Health Clinic   20
## 4           Health Post   86
## 5 Community Health Unit    2
## bar plot of reclassified facility types 
plot_RC_facility_types_WHO = ggplot(RC_facility_types_WHO, aes(x=Var1, y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_RC_facility_types_WHO = plot_RC_facility_types_WHO + labs(x = "Reclassified facility types", y = "Frequency") + coord_flip() + theme_minimal() + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold")) + ggtitle("WHO") + expand_limits(y=c(0,600))
plot_RC_facility_types_WHO

# Types of ownership
ownership_WHO = as.data.frame(table(sf_malawi_WHO$Ownership))
ownership_WHO
##              Var1 Freq
## 1             MoH  463
## 2             FBO  168
## 3 Local authority    5
## 4             NGO    3
## bar plot of ownership
plot_ownership_WHO = ggplot(ownership_WHO, aes(x=Var1, y=Freq)) + geom_bar(stat = "identity", fill="slategray")
plot_ownership_WHO = plot_ownership_WHO + labs(x="Ownership", y = "Frequency") + coord_flip() + theme_minimal() + ggtitle("WHO") + theme(axis.title.y = element_blank(), plot.title = element_text(face = "bold")) + expand_limits(y=c(0,600))
plot_ownership_WHO

Both data sources contain no information on services available, capacity or equipment. MFL does state whether facility is functional.

Classification of MFL facilities aligns more with the structure of the health care system in Malawi (community, primary, secondary, tertiary), it differentiates central hospitals from district and other hospitals. WHO has additional rural and mission hospitals, where do they fit in?

https://www.health.gov.mw/index.php/2016-01-06-19-58-23/national-aids states that at community level, health posts, dispensaries and maternity clinics offer services. Primary includes health centers and community hospitals, secondary consists of district and some CHAM hospitals, tertiary includes central hospitals.

Analysis:

## tmap mode set to interactive viewing
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data
## Warning in sf::st_is_longlat(shp2): bounding box has potentially an invalid
## value range for longlat data

Qs to address?:

  1. How many facilities are in the same location across the MFL and WHO datasets?
  1. Do same facilities share same names/other attributes?
# Qs 1 - how many intersect?

# st_intersection
intersect_WHO_MFL = st_intersection(x=sf_malawi_WHO, y=malawi_facilities_MFL)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
intersect_WHO_MFL ## only 2 intersect directly, so are same up to 5 decimal places?
## Simple feature collection with 2 features and 17 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 33.29456 ymin: -11.53894 xmax: 33.41925 ymax: -11.45836
## geographic CRS: WGS 84
## # A tibble: 2 x 18
##   Country Admin1 Facility.name Facility.type Ownership LL.source iso3c
## * <chr>   <chr>  <chr>         <fct>         <fct>     <chr>     <chr>
## 1 Malawi  North… Euthini Heal… Health Centre MoH       GPS       MWI  
## 2 Malawi  North… Madede Healt… Health Centre MoH       GPS       MWI  
## # … with 11 more variables: facility_type_9 <fct>, CODE <chr>, NAME <chr>,
## #   COMMON.NAME <chr>, OWNERSHIP <fct>, TYPE <fct>, STATUS <chr>, ZONE <chr>,
## #   DISTRICT <chr>, DATE.OPENED <chr>, geometry <POINT [°]>
  1. The ones that aren’t, are they within 50m of another facility?
# Try merge_points()
merge_points("malawi", datasources=c('who', 'healthsites'), dist_same_m=50)
## although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
## malawi num points: who : 639 healthsites : 266  shared at dist thresh 50 m : 52  after merging: 853
  1. How many facilities per admin1 and admin2 regions?
# Facilities per region 
# MFL, admin1
facility_admin1_MFL = st_intersects(malawi_admin1, malawi_facilities_MFL, sparse = TRUE)
malawi_admin1$facility = lengths(facility_admin1_MFL)
head(malawi_admin1) # gave no. of facilities per region
## Simple feature collection with 3 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 32.67036 ymin: -17.12952 xmax: 35.91857 ymax: -9.368326
## geographic CRS: WGS 84
##         shapeName shapeISO           shapeID shapeGroup shapeType
## 1  Central Region     MW-C MWI-ADM1-3_0_0-B1        MWI      ADM1
## 2 Southern Region     MW-S MWI-ADM1-3_0_0-B2        MWI      ADM1
## 3 Northern Region     MW-N MWI-ADM1-3_0_0-B3        MWI      ADM1
##                         geometry facility
## 1 MULTIPOLYGON (((33.5486 -12...      508
## 2 MULTIPOLYGON (((34.58346 -1...      609
## 3 MULTIPOLYGON (((34.02092 -1...      248
# Intersect opposite to get admin1 region for each facility 
admin1_facility_MFL = unlist(st_intersects(malawi_facilities_MFL, malawi_admin1, sparse = TRUE))

# error when adding admin1_facility_MFL data to malawi_facilities_MFL, 
# no. of rows are less in admin1_, why?


# trying to figure out which points don't intersect 

intersect_admin1_MFL = st_intersection(malawi_facilities_MFL, malawi_admin1) # returns data frame of matches

# working out which facilities did not intersect
intersect_admin1_MFL = as.data.frame(intersect_admin1_MFL) # convert to data frame to use in anti_join

no_intersect_admin1_MFL = anti_join(malawi_facilities_MFL, intersect_admin1_MFL) # shows which did not intersect 

# turn back into sf 
no_intersect_admin1_MFL = st_as_sf(no_intersect_admin1_MFL) 
intersect_admin1_MFL = st_as_sf(intersect_admin1_MFL) 

# !!! something is wrong, facilities with same district show different 
# value for region. Was data for district inputted incorrectly or are 
# coordinates wrong?


## Trying to find the source of the error 

## adding admin1_facility_MFL was the issue (the results from st_intersect),
## st_intersect and st_intersection are not computed in the same order
## just use intersect_admin1_MFL for the plots, don't need to combine with facilities that did not
## intersect


# Map with admin1 layer 
tmap_admin1_MFL = tmap_facilities_MFL + tm_shape(st_geometry(malawi_admin1)) + tm_borders()
tmap_admin1_MFL
## Map - number of facilities per region
tmap_admin1_MFL2 = tm_shape(st_geometry(malawi_admin1)) + tm_borders() + tm_shape(malawi_admin1) + tm_fill("facility", style = "cat", palette = sequential_hcl(4, palette = "YlGnBu", rev = TRUE), title = "MFL facilities")

tmap_admin1_MFL2
# Plot of facilities in each admin1 region 

## transform columns
intersect_admin1_MFL$shapeName = as.factor(intersect_admin1_MFL$shapeName)
intersect_admin1_MFL$TYPE = as.factor(intersect_admin1_MFL$TYPE)

## freq of facility types by region
df_facility_admin1_MFL = as.data.frame(table(intersect_admin1_MFL$shapeName, 
                                             intersect_admin1_MFL$TYPE))

## subsetting by region for plots 

# region 1 = Central 
admin1_region1 = filter(df_facility_admin1_MFL, Var1 == "Central Region")

admin1_region1_plot = ggplot(admin1_region1, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "MFL Central Region (total facilities=1365)") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 250))

# region 2 = Southern
admin1_region2 = filter(df_facility_admin1_MFL, Var1 == "Southern Region")

admin1_region2_plot = ggplot(admin1_region2, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "MFL Southern Region (total facilities=1365)") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0,250))

# region 3 = Northern
admin1_region3 = filter(df_facility_admin1_MFL, Var1 == "Northern Region")

admin1_region3_plot = ggplot(admin1_region3, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "MFL Northern Region (total facilities=1365)") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 150))


# MFL, admin2 
facility_admin2_MFL = st_intersects(malawi_admin2, malawi_facilities_MFL, sparse = TRUE)
malawi_admin2$facility = lengths(facility_admin2_MFL)

sum(malawi_admin2$facility) # also has 1365 intersections
## [1] 1365
# Map - number of facilities per district
tmap_admin2_MFL = tm_shape(st_geometry(malawi_admin2)) + tm_borders() + tm_shape(malawi_admin2) + tm_fill("facility", palette = sequential_hcl(6, palette = "YlGnBu", rev = TRUE), title = "MFL facilities") + tm_layout(title = "MFL")

tmap_admin2_MFL
# district for each facility, intersect the other way
intersect_admin2_MFL = st_intersection(malawi_facilities_MFL, malawi_admin2)


# WHO, admin1 (dataset already has admin1 column)

# Map with admin1 layer
tmap_admin1_WHO = tmap_facilities_WHO + tm_shape(st_geometry(malawi_admin1)) + tm_borders()
tmap_admin1_WHO
## intersecting for Map
facility_admin1_WHO = st_intersects(malawi_admin1, sf_malawi_WHO, sparse = TRUE)
malawi_admin1$facility_WHO = lengths(facility_admin1_WHO)
head(malawi_admin1)
## Simple feature collection with 3 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 32.67036 ymin: -17.12952 xmax: 35.91857 ymax: -9.368326
## geographic CRS: WGS 84
##         shapeName shapeISO           shapeID shapeGroup shapeType
## 1  Central Region     MW-C MWI-ADM1-3_0_0-B1        MWI      ADM1
## 2 Southern Region     MW-S MWI-ADM1-3_0_0-B2        MWI      ADM1
## 3 Northern Region     MW-N MWI-ADM1-3_0_0-B3        MWI      ADM1
##                         geometry facility facility_WHO
## 1 MULTIPOLYGON (((33.5486 -12...      508          238
## 2 MULTIPOLYGON (((34.58346 -1...      609          271
## 3 MULTIPOLYGON (((34.02092 -1...      248          130
## Map - number of facilities per region
tmap_admin1_WHO2 = tm_shape(st_geometry(malawi_admin1)) + tm_borders() + tm_shape(malawi_admin1) + tm_fill("facility_WHO", style = "cat", title = "WHO facilities", palette = sequential_hcl(4, palette = "YlGnBu", rev = TRUE))

tmap_admin1_WHO2
# freq of facility types by region
malawi_WHO$Admin1 = as.factor(malawi_WHO$Admin1)
df_facility_admin1_WHO = as.data.frame(table(malawi_WHO$Admin1, 
                                             malawi_WHO$`Facility type`))

# Central
admin1_central = filter(df_facility_admin1_WHO, Var1 == "Central")

admin1_central_plot = ggplot(admin1_central, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "WHO Central Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 250))

# Southern
admin1_southern = filter(df_facility_admin1_WHO, Var1 == "Southern")

admin1_southern_plot = ggplot(admin1_southern, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "WHO Southern Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0,250))

# Northern
admin1_northern = filter(df_facility_admin1_WHO, Var1 == "Northern")

admin1_northern_plot = ggplot(admin1_northern, aes(x=Var2, y=Freq)) + geom_bar(stat = "identity", fill="slategray") + labs(x="Facility type", y="Frequency", title = "WHO Northern Region") + theme(plot.title = element_text(face = "bold")) + theme_minimal() + coord_flip() + expand_limits(y=c(0, 150))


# WHO, admin2 
facility_admin2_WHO = st_intersects(malawi_admin2, sf_malawi_WHO, sparse = TRUE)
malawi_admin2$facility_WHO = lengths(facility_admin2_WHO)

# Map - number of facilities per district 
tmap_admin2_WHO = tm_shape(st_geometry(malawi_admin2)) + tm_borders() + tm_shape(malawi_admin2) + tm_fill("facility_WHO", palette = sequential_hcl(4, palette = "YlGnBu", rev = TRUE), title = "WHO facilities")

tmap_admin2_WHO
# district for each facility, intersect the other way 
intersect_admin2_WHO = st_intersection(sf_malawi_WHO, malawi_admin2)


# final plots
plot_admin1_central = align_plots(admin1_region1_plot, admin1_central_plot, align = "v")
ggdraw(plot_admin1_central[[1]])

ggdraw(plot_admin1_central[[2]])

plot_admin1_southern = align_plots(admin1_region2_plot, admin1_southern_plot, align = "v")
ggdraw(plot_admin1_southern[[1]])

ggdraw(plot_admin1_southern[[2]])

plot_admin1_northern = align_plots(admin1_region3_plot, admin1_northern_plot, align = "v")
ggdraw(plot_admin1_northern[[1]])

ggdraw(plot_admin1_northern[[2]])

  1. Proportion of facility types per district
# district per facility (does not include the 62 facilities) MFL
intersect_admin2_MFL = st_intersection(malawi_facilities_MFL, malawi_admin2)
## Warning in st_is_longlat(x): bounding box has potentially an invalid value range
## for longlat data
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
intersect_admin2_MFL$TYPE = as.factor(intersect_admin2_MFL$TYPE)
intersect_admin2_MFL$shapeName = as.factor(intersect_admin2_MFL$shapeName)

# types by district
by_district_MFL = as.data.frame(table(intersect_admin2_MFL$TYPE, intersect_admin2_MFL$shapeName))
head(by_district_MFL)
##                Var1   Var2 Freq
## 1  Central Hospital Balaka    0
## 2 District Hospital Balaka    1
## 3          Hospital Balaka    0
## 4     Health Centre Balaka   11
## 5            Clinic Balaka   10
## 6       Health Post Balaka    4
# WHO 
intersect_admin2_WHO = st_intersection(sf_malawi_WHO, malawi_admin2)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
intersect_admin2_WHO$Facility.type = as.factor(intersect_admin2_WHO$Facility.type)
intersect_admin2_WHO$shapeName = as.factor(intersect_admin2_WHO$shapeName)

by_district_WHO = as.data.frame(table(intersect_admin2_WHO$Facility.type, intersect_admin2_WHO$shapeName))
head(by_district_WHO)
##                 Var1   Var2 Freq
## 1   Central Hospital Balaka    0
## 2  District Hospital Balaka    1
## 3   Mission Hospital Balaka    0
## 4     Rural Hospital Balaka    0
## 5 Community Hospital Balaka    0
## 6      Health Centre Balaka   13
# Example pie charts for district 

## Balaka 
balaka_MFL = filter(by_district_MFL, Var2 == "Balaka")
balaka_MFL = filter(balaka_MFL, Freq > 0)
balaka_MFL = balaka_MFL[ , -c(2)]

## plot
balaka_plot_MFL = sunburst(balaka_MFL, count = TRUE, legend = list(w=120))
balaka_plot_MFL
Legend
## 
balaka_WHO = filter(by_district_WHO, Var2 == "Balaka")
balaka_WHO = filter(balaka_WHO, Freq > 0)
balaka_WHO = balaka_WHO[ , -c(2)]

## plot
balaka_plot_WHO = sund2b(balaka_WHO, rootLabel = "Facilities")
balaka_plot_WHO
## Lilongwe 
lilongwe_MFL = filter(by_district_MFL, Var2 == "Lilongwe")
lilongwe_MFL = filter(lilongwe_MFL, Freq > 0)
lilongwe_MFL = lilongwe_MFL[ , -c(2)]

## plot
lilongwe_plot_MFL = sund2b(lilongwe_MFL, rootLabel = "Facilities")
lilongwe_plot_MFL
## 
lilongwe_WHO = filter(by_district_WHO, Var2 == "Lilongwe")
lilongwe_WHO = filter(lilongwe_WHO, Freq > 0)
lilongwe_WHO = lilongwe_WHO[ , -c(2)]

## plot
lilongwe_plot_WHO = sund2b(lilongwe_WHO, rootLabel = "Facilities")
lilongwe_plot_WHO
  1. Automating the process of sorting facility types into 4 tiers as (Falchetta et al., 2020). Tier 1: dispensary or basic health post; Tier 2: health center; Tier 3: provincial hospital or regional hospital; and Tier 4: central hospital or national hospital.

Ways to apply the Falchetta classification:

# Falchetta classification with Malawi MFL 
for(i in 1:nrow(malawi_facilities_MFL)) {
   if(malawi_facilities_MFL$TYPE[i] == "Central Hospital"){
      malawi_facilities_MFL[i, "Tier"] = 4
   }
   else if(malawi_facilities_MFL$TYPE[i] == "District Hospital" | malawi_facilities_MFL$TYPE[i] == "Hospital"){
      malawi_facilities_MFL[i, "Tier"] = 3
   }
   else if(malawi_facilities_MFL$TYPE[i] == "Health Centre" | malawi_facilities_MFL$TYPE[i] == "Clinic"){
      malawi_facilities_MFL[i, "Tier"] = 2
   }
   else if(malawi_facilities_MFL$TYPE[i] == "Dispensary" | malawi_facilities_MFL$TYPE[i] == "Health Post"){
      malawi_facilities_MFL[i, "Tier"] = 1
   }
   else{
      malawi_facilities_MFL[i, "Tier"] = 0
   }
}

# With Zambia MFL 
zambia_MFL <- read_csv("https://raw.githubusercontent.com/MOH-Zambia/MFL/master/geography/data/facility_list.csv")
zambia_MFL$facility_type = as.factor(zambia_MFL$facility_type)
levels(zambia_MFL$facility_type)
## [1] "Border Health Post"                "Health Post"                      
## [3] "Hospital - Level 1"                "Hospital - Level 2"               
## [5] "Hospital - Level 3"                "Hospital Affiliated Health Centre"
## [7] "Rural Health Centre"               "Urban Health Centre"              
## [9] "Zonal Health Centre"
# https://healthmarketinnovations.org/sites/default/files/Final_%20CHMI%20Zambia%20profile.pdf states 3 broad
# levels,
# Primary = health centres and health posts
# Secondary = Provincial and district hospitals
# Tertiary = Central hospitals
# https://www.severemalaria.org/countries/zambia/zambia-health-system - describes hospital levels

for(i in 1:nrow(zambia_MFL)) {
   if(is.na(zambia_MFL$facility_type[i]) == TRUE){
      zambia_MFL[i, "Tier"] = 0
   }
   else if(zambia_MFL$facility_type[i] == "Hospital - Level 3"){
      zambia_MFL[i, "Tier"] = 4
   }
   else if(zambia_MFL$facility_type[i] == "Hospital - Level 2" | zambia_MFL$facility_type[i] == "Hospital - Level 1"){
      zambia_MFL[i, "Tier"] = 3
   }
   else if(zambia_MFL$facility_type[i] == "Hospital Affiliated Health Centre" | zambia_MFL$facility_type[i] == "Rural Health Centre" |zambia_MFL$facility_type[i] == "Zonal Health Centre" | zambia_MFL$facility_type[i] == "Urban Health Centre"){
      zambia_MFL[i, "Tier"] = 2
   }
   else if(zambia_MFL$facility_type[i] == "Border Health Post" | zambia_MFL$facility_type[i] == "Health Post"){
      zambia_MFL[i, "Tier"] = 1
   }
   else{
      zambia_MFL[i, "Tier"] = 0 
   }
}

# no. of facility types in each tier across SSA
falchetta_tiers = read_excel("~/malawi-health-facilities-1/parser_healthcare_types.xlsx")
falchetta_tiers = falchetta_tiers[ , -1]

## subset
falchetta_tier4 = filter(falchetta_tiers, Tier == 4)
falchetta_tier3 = filter(falchetta_tiers, Tier == 3)
falchetta_tier2 = filter(falchetta_tiers, Tier == 2)
falchetta_tier1 = filter(falchetta_tiers, Tier == 1)

## select unique elements in types column 
falchetta_tier4_types = unique(falchetta_tier4$ft) # 31 type names 
falchetta_tier3_types = unique(falchetta_tier3$ft) # 62 
falchetta_tier2_types = unique(falchetta_tier2$ft) # 57
falchetta_tier1_types = unique(falchetta_tier1$ft) # 31

# Are facility types ever put into a different tier between countries? 

# Tier 4
tier4_types_check = function(tier_types){ 
  
   tier4_in_other_tiers = c(rbind())
   
   for(i in tier_types) {
     if(i %in% falchetta_tier3$ft == TRUE){
     tier4_in_other_tiers = falchetta_tier3[which(i == falchetta_tier3$ft), ]
     }
     else if(i %in% falchetta_tier2$ft == TRUE){
      tier4_in_other_tiers = falchetta_tier2[which(i == falchetta_tier2$ft), ]
     }
     else if(i %in% falchetta_tier1$ft == TRUE){
      tier4_in_other_tiers = falchetta_tier1[which(i == falchetta_tier1$ft), ]
     }
     else{
      print("No tier 4 facility types present in other tiers")
     }
   return(tier4_in_other_tiers)
   }
}

tier4_types_check(falchetta_tier4_types)
## [1] "No tier 4 facility types present in other tiers"
## NULL
# Tier 3
tier3_types_check = function(tier_types){ 
  
   tier3_in_other_tiers = c(rbind())
   
   for(i in tier_types) {
     if(i %in% falchetta_tier4$ft == TRUE){
     tier3_in_other_tiers = falchetta_tier4[which(i == falchetta_tier4$ft), ]
     }
     else if(i %in% falchetta_tier2$ft == TRUE){
      tier3_in_other_tiers = falchetta_tier2[which(i == falchetta_tier2$ft), ]
     }
     else if(i %in% falchetta_tier1$ft == TRUE){
      tier3_in_other_tiers = falchetta_tier1[which(i == falchetta_tier1$ft), ]
     }
     else{
      print("No tier 3 facility types present in other tiers")
     }
   return(tier3_in_other_tiers)
   }
}

tier3_types_check(falchetta_tier3_types)
## [1] "No tier 3 facility types present in other tiers"
## NULL
# Tier 2
tier2_types_check = function(tier_types){ 
  
   tier2_in_other_tiers = c(rbind())
   
   for(i in tier_types) {
     if(i %in% falchetta_tier3$ft == TRUE){
     tier2_in_other_tiers = falchetta_tier3[which(i == falchetta_tier3$ft), ]
     }
     else if(i %in% falchetta_tier4$ft == TRUE){
      tier2_in_other_tiers = falchetta_tier4[which(i == falchetta_tier4$ft), ]
     }
     else if(i %in% falchetta_tier1$ft == TRUE){
      tier2_in_other_tiers = falchetta_tier1[which(i == falchetta_tier1$ft), ]
     }
     else{
      print("No tier 2 facility types present in other tiers")
     }
   return(tier2_in_other_tiers)
   }
}

tier2_types_check(falchetta_tier2_types)
## [1] "No tier 2 facility types present in other tiers"
## NULL
# Tier 1
tier1_types_check = function(tier_types){ 
  
   tier1_in_other_tiers = c(rbind())
   
   for(i in tier_types) {
     if(i %in% falchetta_tier3$ft == TRUE){
     tier1_in_other_tiers = falchetta_tier3[which(i == falchetta_tier3$ft), ]
     }
     else if(i %in% falchetta_tier2$ft == TRUE){
      tier1_in_other_tiers = falchetta_tier2[which(i == falchetta_tier2$ft), ]
     }
     else if(i %in% falchetta_tier4$ft == TRUE){
      tier1_in_other_tiers = falchetta_tier4[which(i == falchetta_tier4$ft), ]
     }
     else{
      print("No tier 1 facility types present in other tiers")
     }
   return(tier1_in_other_tiers)
   }
}

tier1_types_check(falchetta_tier1_types)
## [1] "No tier 1 facility types present in other tiers"
## NULL
# using left-join 

## table I made, similar to parser_healthcare_types, only for Malawi
Facilities = c("Dispensary", "Health Post", "Clinic", "Health Centre", "Hospital", "District Hospital", "Central Hospital")
Country = c(replicate(7, "Malawi"))
Tier = c(1, 1, 2, 2, 3, 3, 4)

Tiers_MFL = data.frame(Facilities, Country, Tier)
Tiers_MFL
##          Facilities Country Tier
## 1        Dispensary  Malawi    1
## 2       Health Post  Malawi    1
## 3            Clinic  Malawi    2
## 4     Health Centre  Malawi    2
## 5          Hospital  Malawi    3
## 6 District Hospital  Malawi    3
## 7  Central Hospital  Malawi    4
## left_join 
join_MFL = left_join(malawi_facilities_MFL, Tiers_MFL, by = c("TYPE" = "Facilities"))